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with  a  good  attitude. 

It  goes  without  saying  that  every  married  person  that 
goes  to  AFIT  puts  their  spouse  through  a  great  deal  of 
heartache  and  pain.  My  wife  has  endured  two  complete  AFIT 
tours  of  18  months  each  and  deserves  to  be  congratulated  for 
being  an  "AFIT  Widow"  twice  and  still  staying  married  to  me. 


Robert  F.  Stierwalt 


Table  of  Contents 


Preface . 

List  of  Figures . 

Abstract  . 

I.  Introduction  . 

Theory  . 

Research  Proposal  . 

II.  Theory  of  Incoherent  Imaging . 

III.  Superresolution  Scheme . 

Math  Model . 

Properties  of  Ill-conditioned  Matrices 

Least  Squares  Solution . 

Modifications  to  Least  Squares  Method. 

IV.  Computer  Model . 

Introduction . 

Computer  Model . 

V.  Results  and  Conclusions  . 

Results  . 

Band-Pass  vs  Low  Pass  Pupil  .  .  . 
Effect  of  A  Priori  Knowledge  .  . 

Effect  of  a  . 

Effect  of  Noise  . 

Error  Analysis . 

Conclusions  . 

Appendix  A:  Computer  Program  Listing  for 

Superresolution  Program  SRES.  .  .  . 

Appendix  B:  Data  from  Computer  Program  SRES.  . 

Appendix  C:  Typical  Data  from  SRES 

Computer  Program  . 


Bibliography 


Figure 

1.  Impulse  Response . 

2.  Minimum  Rayleigh  Resolution  . 

3.  Generalized  Image  Model  . 

4.  OTF  of  P  »  Rect(x/a) . 

5.  Superresolution  Flow  Chart . 

6.  Examples  of  Modulation  Transfer  Functions 

7.  Superresolution  Example  #1 . 

8.  Superresolution  Example  #2 . 

9.  Superresolution  Example  #3 . 

10.  Superresolution  Example  #4 . 

11.  Error  vs  Known  Length  of  Object  . 

12.  Error  vs  Size  of  Pupil . 


AFIT/GEO/ENP/85D-5 


Abstract 

This  thesis  discusses  the  problem  of  incoherent  imaging 
in  a  diffraction  limited  optical  system.  The  purpose  of  the 
thesis  was  to  prove  that  resolving  two  incoherent  point 
sources  of  light  is  possible  and  achievable  under  certain 
circumstances.  The  effects  of  noise  are  considered  when 
trying  to  superresolve  the  two  Incoherent  objects. 

The  analysis  assumes  a  finite  object  of  known  maximum 
extent  with  an  estimation  of  the  noise  In  the  system.  The 
noise  is  assumed  to  be  Gaussian,  white,  and  additive  for  all 
spatial  frequencies.  The  superresolution  process  uses  the 
standard  least  squares  process  to  achieve  minimum  error  with 
a  smoothing  or  regularization  procedure.  The  singular  values 
of  the  transfer  matrix  are  modified  to  attenuate  the  very 
small  singular  values  to  avoid  noise  amplification  in  the 
high  order  terms.  The  effect  of  the  noise  is  overcome  by  the 
use  of  a  smoothing  parameter,  Q  ,  as  shown  in  the  results. 

The  superresolution  process  works  extremely  well  when  the 
extent  of  the  object  is  known  a  priori  to  have  a  certain 
bound  or  maximum.  Components  of  the  restored  or  processed 
object  outside  the  known  bounds  are  attenuated.  The  results 
indicate  that  band-pass  pupils  can  superresolve  with  only 
limited  knowledge  of  the  object  when  the  smoothing  parameter 
is  used. 
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SUPERRESOLUTION  USING  INCOHERENT  LIGHT 
AND  THE  LEAST  SQUARES  METHOD 

I .  Introduction 

The  purpose  of  this  thesis  is  to  demonstrate  that 
superresolution  of  two  incoherent  point  sources  is  possible 
under  certain  conditions.  Superresolution  or  resolving 
objects  beyond  classical  limits  is  a  current  problem  in 
optics.  It  is  generally  thought  that  diffraction  effects 
represent  the  fundamental  limits  to  optical  system 
performance.  Typically,  objects  placed  closer  than  the 
classical  Rayleigh  limit  (5:309)  cannot  be  distinguished  as 
distinct  objects  but  image  as  a  composite  form.  The  physical 
dimensions  of  lens  systems  traditionally  determine  the 
ultimate  resolving  power  of  the  system. 

In  many  scientific  disciplines  (2:496)  such  as 
spectroscopy  and  astronomy,  superresolution  could  enhance 
research  considerably.  In  general,  for  a  finite  object  in  a 
diffraction  limited  imaging  system,  the  inverse  of  the  linear 
imaging  process  can  be  used  to  yield  the  object.  In  one 
specific  example,  an  object  of  finite  extent  with  known 
maximum  dimensions  can  be  superresolved  to  20%  of  the 
traditional  Rayleigh  criterion  for  two  incoherent  point 
sources  in  a  noise  free  system  (3).  However,  in  all  cases 
the  imaging  process  is  noisy  to  some  degree  and  usually 


produces  a  set  of  ill-conditioned  linear  equations,  thus 
making  the  inversion  of  the  imaging  process  suspect  (4:216). 
Initially,  the  background  and  nature  of  incoherent  imaging 
will  be  presented  using  a  linear  systems  approach.  Following 
the  discussion  of  incoherent  imaging,  a  specific  plan  for 
modeling  a  particular  superresolution  scheme  will  be 
presented  as  a  research  effort  in  order  to  determine  if 
superresolution  is  practical  and  achievable  under  normal 
laboratory  conditions. 

Theory 

Classically  speaking,  the  normal  limit  to  resolution  for 
diffraction  limited  images  is  the  Rayleigh  criterion. 

Incoherent  point  sources  image  as  a  sine  squared  function, 
sinc(x)  =  (sin  ttx)  /  (tx),  in  a  one  dimensional  system  (1:62). 
Figure  1  illustrates  the  1-D  image  for  an  object  consisting 


Figure  1 .  Impulse  Response 


of  an  incoherent  point  source.  The  image  in  this  case  is 
known  as  the  impulse  response  of  the  imaging  system.  The 
dashed  line  represents  the  object  and  the  solid  line 
represents  the  diffraction  limited  image.  Also,  for 
incoherent  light,  the  effect  is  additive  in  the  lrradiance 
distribution  (amplitude  squared)  for  two  or  more  point 
sources . 

For  two  adjacent,  incoherent  point  sources,  separated  by 
the  Rayleigh  criteria,  the  irradiance  distribution  will  be  as 
shown  in  Figure  2.  The  solid  line  represents  the  total 
irradiance  distribution  while  the  dashed  lines  represent  the 
individual  Irradiance  distributions.  Rayleigh  defined  the 
limit  of  resolution  for  circular  apertures  as  the  location  of 
the  first  principle  minimum  of  one  irradiance  distribution 
and  the  first  principle  maximum  of  the  other  irradiance 
distribution  at  the  same  point  in  the  image  plane.  The  same 
limit  can  be  applied  to  rectangular  apertures.  In  Figure  2, 
the  irradiance  has  been  normalized  for  clarity. 

In  1964  J.L.  Harris  (3)  used  the  concept  of  analytical 
continuation  and  a  prior  (known  in  advance)  information  about 
the  object  to  extend  resolution  beyond  the  Rayleigh  limit 
(5:309).  Harris  proved  that  contlnueing  or  extrapolating  the 
function  beyond  the  known  bounds  was  possible  using  the  fact 
that  analytic  functions  are  unique  beyond  the  cutoff 
frequency  of  the  filter  used.  He  proved  that  sampling  the 
Irradiance  distribution  from  a  finite  object  of  known  extent 
could  be  used  along  with  the  properties  of  analytic  functions 


(1:133)  to  resolve  two  point  sources  separated  by  20%  of  the 
classical  Rayleigh  resolution  limit.  This  method  of  Harris 
used  spectral  extrapolation  to  reconstruct  the  object.  Since 
the  two  point  sources  were 


Figure  2.  Minimum  Rayleigh  Resolution 


closer  than  the  Rayleigh  criteria,  this  reconstruction  scheme 
was  also  a  superresolution  process.  This  theory  did  not  hold 
up  under  severe  scrutiny  as  Harris  neglected  noise  and  did 
not  know  at  the  time  that  his  set  of  simultaneous  equations 
was  extremely  ill-conditioned  (3:1481). 

After  analytical  continuation  was  used  to  extrapolate 
the  spectral  components  of  a  finite  object  with  a  priori 
knowledge  of  its  maximum  dimensions,  it  was  realized  that 
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noise  must  be  taken  into  account  (6,7)  to  achieve  any 
measurable  degree  of  superresolution.  One  of  the  earliest 
methods  used  to  superresolve  objects  was  the  iteration  method 
proposed  by  Gerchberg  (7).  Gerchberg  utilized  a  specific 
iteration  method  derived  from  the  general  method  of  Youla 
(8).  The  method  of  Youla  uses  orthogonal  projections  in  a 
well-posed  Hilbert  space  and  a  priori  knowledge  of  the 
extent  of  the  object  to  remove  any  spectral  component  higher 
than  the  cut-off  frequency.  Gerchberg  theorized  that  any 
component  measured  at  a  higher  frequency  than  the  cut-off 
frequency  had  to  be  a  noise  component,  and  was  therefore 
subtracted  out.  This  method  was  iterated  until  the  measured 
output  past  the  cut-off  frequency  was  below  some  threshold 
level.  Gerchberg  also  reasoned  that  noise  could  not  be 
analytically  continued  since  it  was  not  band  limited. 

Gerchberg  did  not  realize  at  the  time  that  his  matrix 
methods  were  unstable  as  Byrne,  et  al,  pointed  out  in  their 
1983  article  (6).  Broadband  noise  served  to  cause  spurious 
oscillations  in  the  data  making  the  restoration  scheme 
suspect.  After  a  good  noise  estimation  was  used,  Byrne,  et 
al ,  were  able  to  use  a  Gerchberg  type  algorithm  to  obtain  a 
reasonable  superresolved  object. 

To  further  the  work  of  superresolution,  Mammone  and 
Eichman  (2)  used  optimum  linear  programming  techniques  to 
smooth  the  data  thus  providing  a  stable,  well-conditioned 
matrix  solution  to  the  image  restoration  problem.  Their 
initial  assumption  was  sirailiar  to  that  of  Harris  (3)  in  that 
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the  Image  irradlance  could  be  sampled  and  the  object 
reconstructed  from  a  set  of  linear  equations.  The  matrix 
solution  can  be  represented  as 

~i  =  Ao  (1.1) 

a  A 

where  1  and  o  are  n  dimensional  vectors  representing  the 
image  irradlance  and  object  irradiance  respectively.  A  is  an 
n  by  n  transformation  matrix  obtained  from  the  discretized 
solution  of  the  object  vector  o  from 

i(xi)  =  |h(xi;xo)|2  *  o(xo)  (1.2) 

The  complete  description  of  Incoherent  imaging  is  contained 
in  the  next  chapter.  Mammone  and  Eichman  chose  to  make  the 
transformation  matrix.  A,  stable  by  filtering  methods.  By 
filtering,  high  frequency  noise  is  eliminated  as  well  as 
making  the  Inversion  of  the  matrix  A  possible,  thus  yielding 
a  solution  for  o. 

While  all  of  the  above  methods  use  spectral 
extrapolation  and  symmetric  low  pass  spectral  filters, 

Cathey,  et  al ,  (9)  have  proposed  a  superresolution  concept 
utilizing  the  same  bandwidth  as  the  extrapolation  method,  but 
using  bandpass,  or  multi-aperture  systems  and  interpolation 
instead  of  extrapolation  to  achieve  superresolution.  It  is 
postulated  that  for  each  imaging  situation,  there  exists  a 
potential  "best"  aperture  window  to  superresolve  the  object. 
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Using  similiar  reconstruction  techniques  in  the  frequency 
domain,  the  object  was  consistently  better  resolved  using 
band  pass  techniques  instead  of  low  pass  filtering. 


Research  Proposal 

In  this  thesis  effort,  incoherent  light  will  be  used 
along  with  a  multiple  aperture  system  to  develop  a 
mathematical  model  for  superresolving  an  object.  By  using 
incoherent  light,  the  transfer  function  of  the  pupil,  or 
aperture,  is  not  as  straight  forward  as  the  coherent  case, 
but  the  matrix  method  of  inverting  the  transformation  matrix 
seems  viable.  To  effectively  evaluate  the  optimum  pupil 
function,  that  pupil  function  yielding  the  minimum  error  will 
be  judged  as  the  best  pupil  function  for  a  given  fill  ratio, 
where  the  fill  ratio  is  the  number  of  apertures  divided  by 
the  total  number  of  available  windows  (9:247).  The  premise 
to  be  proved  is  that  for  a  given  fill  ratio,  band  pass 
filtering  instead  of  low  pass  filtering  provides  better 
resolution  using  the  proposed  algorithm. 


Theory  of  Incoherent 


The  purpose  of  this  chapter  is  to  summarize  the  theory 
of  incoherent  imaging.  This  summary  utilizes  a  linear 
systems  approach  similiar  to  chapter  six  of  reference  one. 

By  using  a  systems  approach,  the  properties  of  an  optical 
system  can  be  given  in  terms  of  input,  transfer  function  or 
impulse  response,  and  output,  irrespective  of  the  number  and 
type  of  internal  optical  elements.  For  optical  systems, 
these  properties  can  be  summarized  in  terms  of  its  exit  or 
entrance  pupil,  effective  focal  length,  and  its  output  with 
an  input  consisting  of  a  point  source  on  the  optic  axis.  In 
this  particular  case,  images  are  located  as  predicted  by 
geometric  optics  and  the  system  is  considered  to  be 
!•  diffraction  limited  (1:103). 


FIGURE  3.  Generalized  1-D  Image  Model 
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Referring  to  figure  3,  the  imaging  system  can  be 
represented  as  a  "black  box"  consisting  of  an  entrance  pupil 
and  an  exit  pupil.  The  properties  of  the  entire  system  can 
be  completely  described  by  specifying  the  properties  of  the 
entrance  pupil  or  exit  pupil.  It  is  assumed  that  the  entire 
optical  system  can  be  adequately  described  by  geometric 
optics  and  that  all  diffraction  effects  can  be  associated 
with  the  entrance  or  exit  pupil.  The  entrance  and  exit 
pupils  are  geometric  projections  of  one  another  which  enables 
an  equivalent  analysis  using  either  one  (1:102-103). 

Since  geometric  optics  adequately  describe  the  behavior 
of  light  between  the  entrance  and  exit  pupils,  diffraction 
dominates  the  behavior  of  light  from  the  object  to  the 
entrance  pupil  and  from  the  exit  pupil  to  the  geometric  image 
plane.  Since  diffraction  effects  can  be  associated  with 
either  the  entrance  or  exit  pupil,  all  diffraction  effects 
will  be  associated  with  the  exit  pupil  in  the  context  of  this 
paper.  This  approach  is  acceptable  since  the  exit  and 
entrance  pupil  are  geometric  images  of  each  other  (1:102). 

Using  the  notation  found  in  Goodman's  text  (1)  the  image 
amplitude  distribution  can  be  expressed  as 

oo 

Uj(xi)  *  Jh(xi  ;xo)U0  (xo)dxo  (2.1) 

-oo 

where  U, (xi)  is  the  image  amplitude,  h(xi)  is  the  impulse 
response  or  transfer  function,  and  U„ (xo)  is  the  object 


amplitude  distribution.  The  Image  plane  is  xi ,  the  object 
plane  xo,  and  the  exit  pupil  plane  is  x.  If  the  lens  law  for 
imaging  is  satisfied,  then  h(xi)  is  the  Fraunhofer 
diffraction  pattern  of  the  exit  pupil  centered  at  xi  *  -Mxo  , 
where  M  is  the  magnification  of  the  system.  Notice  that  the 
system  is  inverted  as  represented  by  the  negative  sign 
(1:95).  The  impulse  response  can  be  written  as 

oo 

h(xi;xo)  =  k/p  ( x)exp( ( -2ir j/Xdi )  (xi-Mxo)x)dx  (2.2) 

-ao 

with  the  understanding  that  the  pupil  function,  P(x) ,  is 
either  zero  or  one  depending  on  whether  it  blocks  or  passes 
light  in  that  particular  interval,  dx.  In  equation  (2.2), 
h( xi ; xo)  can  be  described  as  the  Fourier  Transform  (FT)  of 
the  pupil  function  evaluated  at  the  spatial  frequencies  fx  = 
xl/Xdi.  Also,  K  is  a  complex  constant  that  will  be  discussed 
later.  By  approplate  change  of  variables,  the  impulse 
response  can  be  rewritten  as 

00 

h' (xi ;XO)  =  (Xdix1  )exp(-27rj(xi-xo'  ))dx'  (2.3) 

-00 

where  x'  =  x/Xdi  and  xo ' =  Mxo.  By  defining  0g(xi)  as  the 
geometric  image,  the  real  image  can  be  described  as  the 
convolution  of  the  geometric  image  with  h'(xi),  the  modified 
impulse  response.  Equation  (2.4)  is  the  convolution  integral 
representing  the  final  image  as  the  convolution  of  the 
geometric  image  with  the  modified  impulse  response,  h'(xi). 


10 


Since  the  final  expression  will  be  normalized  with  respect  to 
magnitude,  all  constant  multipliers,  real  or  complex  can  be 
ignored  ( 1 : 105 ) . 

Up  to  this  point,  amplitude  has  been  the  subject  of  our 
derivations,  but  irradiance  (watts/area)  is  the  measured 
quantity.  The  irradiance  is  the  infinite  time  average  of  the 
amplitude  squared  as  shown  in  equation  (2.5). 

Ii  =  ^Ui  (xi)Uf  (xif>  (2.5) 

The  brackets  denote  the  infinite  time  average  and  Uf(xi) 
designates  the  complex  conjugate  of  U;(xi).  For  real  sources 
and  incoherent  light, 

Ig(xo)  =^U9(xo)^>2  (2.6) 

where  g  designates  geometrically  predicted  quantities.  The 
final  irradiance  image  is 

so 

Ii(xi)  =  ky"|h(xi-xo  1  )|I9  (xo  1  )dxo  1  (2.7) 

-oo 

where  the  k  represents  a  constant  multiplier.  It  can  be 
easily  seen  that  the  final  irradiance  image  is  again  a 
convolution  of  a  transfer  function  and  a  predicted  geometric 


quantity.  The  irradiance  transfer  function  is  the  modulus  of 
the  impulse  response,  h'(xi),  squared.  This  implies  that 
there  is  no  phase  information  for  incoherent  imaging .( 1 : 109 ) 
For  ease  and  convenience,  let 

Gg(fx)  =  FT ( la ( XQ 1 ) )  (2.8) 

FT(Ia(xo')),  evaluated  at  fx=0 

where  Gg(fx)  represents  the  normalized  Fourier  transform  of 
the  geometrically  predicted  irradiance  image  of  the  object. 
Likewise 


i 


Gi  (fx)  =  FT  ( Ii  (Xi)  )  (2.9) 

FT( Ii  (xi ) ) ,  evaluated  at  fx=0 

where  G  (fx)  is  the  normalized  Fourier  transform  of  the 
diffraction  limited  image.  The  irradiance  transfer  function, 
H( fx) ,  is  (1:114) 

H  (  fx )  =  FT(  I  h 1  (xl)l2  )  (2.10) 

FT(|h'(xi)|2  ),  evaluated  at  fx=0 

where  again  normalization  has  taken  place.  Thus,  in  the 
frequency  domain, 


Gi  (fx)  =  H(fx)Gg (fx) 


(2.11) 


This  last  equation  can  be  derived  using  the  convolution 
theorem  (13:314) . 

The  optical  transfer  function,  H,  of  the  incoherent 
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irradiance  imaging  system  is 

H  =  !  FT ( h 1  )|2 


(2.12) 


where 


h'  =  FT ( P )  (2.13) 

It  can  be  shown  that  for  f(x)  a  real  function  of  x  that 

FT(  f  (x)  **  f  ( x )  )  =  |  FT  (  f  ( x)  )|2  (2.14) 

By  using  the  autocorrelation  theorem  (13:200).  Therefore 

|  FT  ( h 1  )  |2  =  P»*P  (2.15) 

which  leads  to  the  final  expression  for  the  optical  transfer 
function,  H: 

H  =  P**P  (2.16) 

where  **  denotes  autocorrelation.  Also,  all  subscripts,  etc, 
have  been  dropped  for  clarity.  Now,  equation  (2.11)  can  be 
expressed  as 


Gi  *  (P*P)Gb 


(2.17) 


which  i.~  the  final  spectral  domain  representation  of  the 
imaging  relationship  between  the  object  and  the  optical 
system.  Figure  4  illustrates  the  normalized  transfer 
function  for  a  rectangular  pupil  function.  As  can  be  seen, 

the  OTF ,  H,  has  a  definite  cutoff  frequency,  fc,  which  is  a 
function  of  the  system  parameters.  The  modulus  of  H,  |  H  | , is 
known  as  the  modulation  transfer  function,  MTF  (1:114). 

For  optical  systems  it  would  seem  that  a  simple  Fourier 
inverse  of  the  frequency  domain  representation  of  the  image 
would  easily  yield  the  original  object.  Equation  (2.18) 
represents  the  inverse  Fourier  transform  needed  to  recover 


the  object  from  the  spectral  parameters 


FT“  1  [Gi] 


(2.18) 


The  problem  with  resolution  is  in  the  use  of  limited 
apertures  and  noise,  since  limited  apertures  attenuate  the 
high  spatial  frequency  terras  essential  for  high  resolution 
and  the  noise  terms  are  greatly  magnified  in  the  inverse 
process,  as  will  be  explained  later.  Larger  pupils  would 
alleviate  much  of  the  problem  but  there  is  a  physical  limit 
to  the  size  of  usable  optical  systems,  especially  for  space 
applications . 

In  the  next  chapter,  a  particular  matrix  method  for 
increasing  resolution  with  limited  spatial  bandwidth  will  be 
examined.  It  will  be  shown  that  increased  resolution  can  be 
achieved  under  certain  conditions.  By  knowing  a  priori  that 
the  object  is  finite  with  a  known  maximum  extent,  and  by 
using  a  smoothing  parameter,  the  increase  in  resolution  can 
be  quite  substantial. 


Introduction 


The  purpose  of  this  chapter  is  to  develop  a  mathematical 
model  for  incoherent  imaging  with  limited  bandwidth  optical 
systems.  The  limited  bandwidth  is  realized  by  symmetric 
multiple  apertures  as  well  as  the  traditional  low  pass  (in 
spatial  frequency)  model.  It  has  been  speculated  that  for  a 
limited  bandwidth  system,  better  resolution  could  be  obtained 
using  a  bandpass  aperture  insead  of  the  low  pass  system 
traditionally  used  (9).  Multiple  aperture  systems  consisting 
of  relatively  small,  precise  optical  elements  could  simulate 
larger  optical  systems  which  are  costly  and  extremely  hard  to 
manufacture.  The  optical  systems  will  be  described  in  terms 
of  its  optical  transfer  function,  H,  and  its  exit  pupil,  P. 
The  mathematical  model  uses  discrete  values  with  vector  and 
matrix  analysis.  For  the  purpose  of  this  thesis,  the 
truncation  errors  and  sampling  errors  will  not  be  discussed. 
The  data  is  assumed  to  be  adequately  sampled  and  the 
truncation  errors  are  considered  to  be  negligible  compared  to 
the  noise.  This  chapter  describes  the  discrete  solution  to 
equation  (2.18)  using  linear  methods. 

Math  Model 

To  effectively  model  the  Imaging  system,  a  priori 

I 

knowledge  of  the  object  is  of  prime  importance  in 
reconstrucing  the  object  from  image  data.  The  object  is 


known  to  be  finite  within  the  first  M  units  of  the  N  bit 


object  vector,  x0 .  The  finite  object  can  be  represented  by 
Dxt  where  6  is  the  diagonal  matrix  (all  elements  zero  not  on 
the  diagonal) 

D  =  dlag  ( 1 , 1 . 1,0,0 . 0 )  (3.1) 

M  N-M 

The  matrix  D  is  known  as  the  spatial  truncation  matrix.  The 
spatial  truncation  matrix  has  been  zero-filled  to  an  N  x  N 
matrix  to  match  the  order  of  subsequent  matrices. 

The  N  bit  discrete  Fourier  transform  (DFT)  of  the 
truncated  object  is 

DFT(xt)  =  FDx  (3.2) 

and  F  is  the  matrix  representing  the  Fourier  transformation, 
DFT,  whose  components  are 

F(m.n)  =  exp(  -  j27rmn/N)  (3.3) 

for  m,n  =  0,N-1 

Thus,  F  is  a  complex,  N  x  N  matrix. 

In  any  optical  optical  system  the  pupil  is  finite  and 
passes  only  a  limited  number  of  spatial  frequency  terms.  In 
an  incoherent  imaging  system,  the  transfer  function  is  the 
autocorrelation  of  the  pupil  function,  as  derived  in  chapter 
two.  Since  only  a  finite  number  of  spatial  frequency  terms 


will  pass  through  the  optical  system,  the  filtering  effect, 
or  transfer  function,  can  be  represented  by  the  diagonal 
matrix  B,  whose  diagonal  components  contain  the  appropiate 
attenuation  factors.  The  matrix  B  is  equal  to 

B  =  diag(DC,f  |  ,  f2 . fc  . 0 , 0 .  .  .  . 0 ,  fc  ,  .  .  ,  f ,  ,  f ,  )  (3.4) 


where  DC  is  the  DC  coefficient,  the  fj  's  are  the  various 
spatial  frequency  coefficients,  and  fc  is  the  cutoff 
frequency  (highest  spatial  frequency  passed  by  the  system) . 
The  symmetry  of  the  B  matrix  is  the  same  as  the  symmetry  of 
the  forthcoming  DFT 1 s .  For  appropiate  multiplication  of 
vectors  and  matrices,  either  the  B  matrix  is  changed  to  the 
symmetry  of  the  DFT  (which  it  is)  or  the  DFT  is  rearranged  to 
the  symmetry  of  th  B  matrix.  Thus,  the  spatial  frequency 
representation  of  the  image  is  BFDXo •  The  diagonal  elements 
of  B  are  obtained  from  the  specific  values  obtained  from 

-k 

calculating  P**P.  In  this  thesis,  all  pupils  will  be 
centered  and  symmetric  about  the  optic  axis.  It  is  a 
property  of  discrete  sequences  that  an  N-bit  sequence 
convolved  with  an  M-bit  sequence  generates  a  sequence  of 
length  M+N-l  (15:12).  Thus,  for  a  pupil  of  length  L,  the 
transfer  function  response  will  be  of  length  2L-1.  As  an 
example,  let  P=(llli)'  (Pisa  column  vector,  so  1 
denotes  transpose  of  the  row  vector).  The  convolution  of  P 
with  itself  will  yield  the  diagonal  elements  of  the  matrix  B, 
since  convolution  with  itself  is  autocorrelation. 
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(3.5) 


where  P**P  has  been  rearranged  to  the  same  symmetry  as  the 
DFT  in  equation  (3.2).  In  the  example,  the  transfer  function 
consists  of  seven  elements  which  correspond  to  the  DC  term 
plus  three  positive  and  three  negative  spatial  frequency 
terms.  Therefore,  the  four  bit  pupil  passes  the  first  three 
spatial  frequency  terms  plus  the  zero  frequency  or  DC  term. 

So  far,  the  image  is  represented  as  BFDjC,  .  The  last  step 
remaining  is  to  convert  the  spectral  representation  of  the 
image  back  to  a  conventional  spatial  representation.  The 
normal,  noise  free,  space  domain  image  can  be  written  as 


X;  =  FT  BFDXc 


(3.7) 


which  physically  can  be  expressed  as 


Xi  =  Ax0  +  n 


(3.8) 


9 


where  A  =  FT''  BFD  and  n  is  the  Gaussian,  white,  and  additive 
noise  vector.  The  components  of  FT-1  (inverse  Fourier 
transform  matrix)  are 

FT~’(m»n)  =  exp(2jrmn/N)  (3.9) 

for  m.n  =  0,N-1 

FT-1  is  also  an  N  x  N  complex  matrix. 

Solving  for  x0  in  equation  (3.8)  is  the  mathematical 

problem  of  the  superresolution  process,  and  thus  is  the  heart 

of  this  thesis  effort.  Equation  (3.8)  represents  the 

linear  transformation  of  the  irradiance  from  the  object  plane 

to  the  image  plane.  Also,  the  problem  is  compounded  by  the 

—a 

fact  that  the  image,  Xj  ,  is  affected  by  noise,  n,  in  the 
imaging  system.  For  this  paper,  the  noise  is  considered  to 
be  Gaussian,  white,  and  additive  for  all  frequencies 
considered.  The  solution  of  x0  from  equation  (3.8)  is  not  a 
trivial  matter  as  the  matrix  A  is  severely  ill-conditioned 
which  can  lead  to  serious  problems  in  the  solution  of  x,.  A 
discussion  of  ill-conditioned  matrices  follows. 

Properties  of  Ill-conditioned  Matrices 

The  matrix  A  can  be  severely  11-conditioned. 
Ill-conditioned  matrices  have  the  potential  to  cause  very 
large  errors  when  used  to  solve  linear  equations  because  of 
the  propogation  of  errors.  The  length,  or  size,  of  a  matrix, 
can  be  expressed  as  its  norm  or  magnitude,  and  is  expressed 
as 
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Two  norm  of  A 


(3.10) 


II  A  ! 

where  the  two  norm  of  A  Is  the  most  restrictive  norm  (11:26). 
The  two  norm  of  a  matrix  is  defined  as  the  largest  singular 
value  (the  singular  values  are  represented  as  3;,  for  i  =  l,k 
where  k  is  the  rank  of  the  matrix)  of  that  matrix.  The 
singular  values  of  a  matrix  are  the  square  roots  of  the 
corresponding  eigenvalues  of  the  matrix.  The  condition 
number  of  a  matrix,  c(A) ,  is  a  measure  of  the  stability  of 
the  matrix.  The  condition  number  of  the  matrix  A  is  defined 
as 


c  ( A)  »  ||A||2  x  1 1  A~'|  1 2  (3.11) 

where  the  magnitude  of  A"1  is  equal  to  the  inverse  of  the 
smallest  singular  value  of  A  (12:166). 

As  an  example,  let  A  be  the  32  x  32  matrix  obtained  from 
(see  Appendix  B) : 

P  =  (1111100000011111)'  and 
D  =  Diag ( 1 ,  . . . 1 ) ,  for  M  =  32 

_ i 

where  '  denotes  transpose  since  P  is  a  column  vector.  From 
the  data  table,  the  largest  singular  value  of  A  is  320  and 
the  smallest  singular  value  is  1.46  E-4.  Since  the  magnitude 
or  two  norm  of  A  is  equal  to  1/(1. 46  E-4),  the  condition 
number  is  roughly  2E+10.  Why  is  the  condition  number 
Important?  The  importance  of  the  condition  number  is 


The  condition  number  is  a  measure  of  stability  m  that 
it  gives  an  upper  bound  on  the  possible  magnitude  of  error  in 

— _k  — *  — 

the  solution  of  the  matrix  equation  Ay  =  b,  for  A  the  image 
matrix,  y  the  object,  and  b  the  noisy  object.  Let  e(y) 
represent  the  potential  error  possible  from  the  solution  of 
y=A’,b,  and  let  e(b)  represent  the  error  present  in  b,  the 
Gaussian  noise.  The  relative  error  in  the  solution  for  y 
will  have  the  upper  bound  expressed  by 


e(y)  <=  c ( A )  x  e( b) 


(3.12) 


For  the  example  cited  above,  the  condition  number  was 
approximately  2E+7 ,  so  the  relative  error  limits  can  be 
expressed  as 


e(y )  <=  2E+7  x  e(b) 


(3.13) 


where  e(b)  is  usually  expressed  as  the  noise  variance.  It  is 
important  to  realize  that  the  limit  to  the  error  is  an  upper 
bound  on  the  error  and  that  not  every  component  of  the 
solution  vector  will  be  in  error  by  this  amount,  but  each 
component  could  be  in  error  by  this  amount.  The  use  of 
condition  numbers  to  test  matrices  for  stability  is  a  figure 
of  merit  type  relationship.  For  good  linear  systems,  the 
condition  number  should  be  small.  By  using  the  two  norm,  the 
most  restrictive  error  potential  was  achieved.  Other, 
simpler  norms  could  have  been  used,  but  generally  lead  to 


higher  condition  numbers.  Another  way  of  looking  at 
condition  numbers  is  to  say  that  a  matrix  with  a  large 
condition  number  is  an  almost  singular  matrix,  meaning  that 
the  matrix  inversion  is  quite  suspect.  From  linear  algebra, 
a  unique  solution  to  a  matrix  transformation  is  only  possible 
when  the  inverse  to  the  matrix  exists.  Therefore,  solving  a 
linear  system  of  equations  with  an  almost  singular  matrix 
means  that  the  matrix  is  not  very  stable  and  can  cause  large 
errors  in  the  solution. 

Least  Squares  Solution 

The  problem  of  imaging  can  be  restated  as 


where  =  implies  that  a  least  squares  solution  is  being 
sought.  The  notation  uses  Xjj  as  the  ith  component  of  the 
vector  Xj  .  The  least  squares  solution  minimizes  the 
difference  expressed  by 

1 1  Ax0  —  x ,•  j  |  2  =  error  (3.15) 

where  the  subscript  2  represents  the  2  norm  of  the  vector, 
also  known  as  the  Euclidean  distance,  as  shown  in  equation 


To  implement  the  least  squares  solution,  the  matrix  A 
must  be  decomposed  into  a  product  of  a  diagonal  matrix  and 
two  orthogonal  matrices.  The  components  of  the  diagonal 
matrix  will  be  the  singular  values  of  A.  The  matrix  A  ,  where 
A  =  FT-'  BFD ,  can  be  rewritten  as  a  product  of  three  matrices 
( 12 : 237) 

A  =  USV'  (3.17) 

where  U  is  an  N  X  N  orthogonal  matrix,  and  V'  is  an  L  X  L 
orthogonal  matrix.  S  is  an  N  X  L  orthogonal  matrix  whose 
diagonal  elements  are  the  singular  values  of  A  and  are  in 
decreasing  order.  S  looks  like 

S  ■  Diag(s,  ,s2,s3 . sm  ,  0 . 0 - 0)  (3.18) 

where  s,  >  s  ,  >  s^>...>sm.  There  are  only  M  singular  values 
since  D  reduces  the  rank  of  A  to  M. 

The  least  squares  method  of  solving  linear  equations 
serves  to  minimize  the  two  norm  of  the  difference  vector  as 
expressed  in  equation  (3.16).  Substitute  for  A  in  (3.15)  to 
obtain 

IIusV'jTo  -  "x i  i  1 2  =  error  (3.19) 

where  A=USV'.  Since  U  is  orthogonal,  UU'  =  I,  where  I  is  the 


identity  matrix.  By  multiplying  each  term  by  U',  the 
expression  becomes 


I  SV'Xo 


error 


(3.20) 


where  the  error  is  the  same,  since  multiplying  a  vector  by  an 
orthogonal  matrix  does  not  magnify  the  norm  of  the 
vector(  12-.282)  .  By  substitution,  equation  (3.20)  becomes 


[  |  Sy  -  b  *1 1  2  =  error  (3.21) 

“  “  — k  k 

where  y=V'x0  and  D'=Uxj .  Therefore,  x0  solves  the  least 
squares  problem  (3.15)  if  and  only  if  y=Vx0 solves 


i  Sy  —  b '  [  1 2  -  minimum  error 


(3.22) 


Since  Sy  is  a  diagonal  matrix  composed  of  the  singular  values 
of  A  times  y,  ,  the  difference  vector  of  (3.15)  can  be 
expressed  as 

*>'  1 12  “  [is,  y,  -  b,  '  )2  +  ...  { smym  -  b,;  )a  +  b^, 

(3.23) 


The  minimum  error  solution  satisfies 


for  all  values  of  i .  By  back  substitution. 


Ym  =  b^  /9m  =  £u'(m,j)  xu  (1/Sj) 

j=i 


and  since  y=V'x0 


(3.25) 


V.V.-l 
.v.v.v 1 
•--.■-■'.'I 


xo  =  Vy 


(3.26) 


since  V  and  V'  are  orthogonal  matrices.  By  substituting  for 
y  in  equation  (3.26),  the  mth  component  of  x«  is 


L  _  L  _ 

(Xo)m  =  Y,  V(m,k)£u'  (m,  j)  (x,j)  (1/Sj) 

K-’l  j:l 


which  can  be  rearranged  to  yield 


(3.27) 


i.  L  _  _ 

(Xo)rn  U  •  (m,  j  )  (x,j)  ( 1/Sj  )V(m,k) 

K-l  J  =  1 


(3.28) 


By  using  vector  notation,  equation  (3.27)  and  (3.28)  can 


be  rewritten  as 


Xo  =  £(1/Sj  )  (x  j  .  uJ  )  vJ 

j  =) 


(3.29) 


where  uJ  and  vd  are  the  left  and  right  singular  vectors 

mm 

obtained  from  the  jth  column  of  the  matrices  U  and  V 
respectively.  The  solution  vector,  x0 ,  denotes  a  least 
squared  solution,  and  •  is  the  dot  product  operator.  This 
solution  is  also  referred  to  as  the  inverse  filtering 
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solution,  since  the  least  squares  solution  serves  to  reverse 
the  matrix  operation  and  return  the  original  (object)  vector 
in  our  specific  application  (14:205). 


Modifications  to  Least  Squares  Method 

Since  there  are  very  small  singular  values  in  the  matrix 
A,  the  presence  of  noise  will  greatly  reduce  the 
effectiveness  of  the  least  squares  solution  in  the  high  order 
terms.  By  using  the  method  of  Rushforth,  et  al,  (14),  l/sK  is 
replaced  with  an  appropiate  smoothing  or  regularization 
function  to  attenuate  the  high  order  terms  and  pass  the  lower 
order  terras  unchanged.  The  smoothing  function  serves  to 
avoid  overaraplif ication  of  the  higher  order  noise  terras.  For 
this  thesis,  l/sK  is  replaced  with  f(sK),  where 

f(sK)  = _ s* _  (3.30) 

SK  +  « 

Choosing  an  appropiate  <*  prevents  overamplification  of  the 
high  order  noise  terms.  “  was  chosen  in  most  cases  to  be 
10E-10  as  this  provided  adequate  attenuation  with  respect  to 
the  smaller  singular  values.  Various  other  smoothing 
functions  and  concepts  could  be  considered.  The  basic 
concept  is  to  reject  those  higher  order  terms  knowing  as  much 
about  the  object  as  possible  beforehand.  To  get  high 
resolution,  some  high  order  terms  must  be  present  in  the 
solution.  Therefore,  if  the  object  is  known  beforehand  to  be 
less  than  10  units  long,  the  spatial  truncation  matrix,  D, 


-r-  1 .  w .  P.  W.'W.'  >.  yjw.n«jy.fy 


solution.  Therefore,  if  the  object  is  known  beforehand  to  be 
less  than  10  units  long,  the  spatial  truncation  matrix,  6, 
will  be  finite  for  only  9  elements.  By  truncating  at  nine 
elements,  the  matrix  A  will  be  of  rank  less  than  or  equal  to 
nine,  so  there  will  be  only  nine  singular  values,  thus  the 
lower  singular  values  will  already  be  rejected  (12:336). 

The  following  chapter  will  describe  the  computer  program 
and  computational  algorithms  used  to  implement  the  least 
squares  solution. 


28 


IV.  Computer  Model 

Introduction 

The  purpose  of  this  chapter  is  to  describe  the 
analytical  implementation  of  the  least  squares 
superresoslution  scheme  as  described  in  the  previous 
chapters.  The  computer  system  used  was  a  VAX  11/785 
utilizing  the  VMS  4.1  operating  system.  The  subroutines  used 
for  convolution,  Gaussian  random  number  generation,  and 
singular  value  decomposition  of  matrices  were  from  the 
International  Mathematical  and  Statistical  Subroutine  Library 
(IMSL)  (10).  The  matrix  multiplication  routines  were  written 
by  the  author.  The  program  was  written  using  standard 
Fortran  77  ( 18)  . 

Figure  5.  is  a  flowchart  depicting  the  superresolution 
process  as  described  in  the  previous  chapter.  Refer  to  this 
figure  throughout  the  chapter  for  reference.  A  complete 
listing  of  the  program,  SRES,  is  contained  in  the  appendix. 


Computer  Model 

The  following  variables  are  defined  for  convenience. 

L  =  pupil  length  (integer) 

N  =  dimension  of  square  matrix  A  (integer) 

M  =  finite,  known  extent  of  object 

=  smoothing  coefficient  (real  exponential) 

P  =  pupil  vector 

DSEED  =  double  precision  constant 
sj  =  Jth  singular  value  (real) 

±  =  convolution  operator 
Q.  =  matrix  named  Q 
r  =  vector  named  r 
r,  =  ith  component  of  r 


A-1  =  inverse  of  matrix  A  _ 

A'  =  transpose  of  matrix  A 
SNR  =  power  signal  to  noise  ratio 
PSEQ  =  power  in  discrete  sequence 
On  =  noise  variance 
•  =  dot  product  operator 

The  first  variable  to  define  is  N,  dimension  of  the  transfer 
matrix,  A.  The  criterion  for  N  was  that  it  be  large  enough 
to  model  an  imaging  system  but  small  enough  to  operate 
efficiently  with  the  computer.  Also,  N  was  chosen  to  be  a 
power  of  two,  since  this  makes  the  system  easier  to  work 
with,  with  respect  to  the  IMSL  subroutines.  N  was  chosen  to 
be  32  for  all  configurations.  By  choosing  a  value  for  N,  the 
values  for  M  and  L  are  somewhat  restricted.  The  length  of 

the  pupil,  L,  must  be  less  than  or  equal  to  N/2,  since  the 

—  — * 

matrix  B  is  composed  of  elements  of  P*P  and  convolution  of 

two  vectors  lenghens  the  resulting  vector.  Also,  the  object 

dimension,  M,  cannot  be  larger  than  N  because  of  the  matrix 

vector  multiplication  scheme. 

The  computer  program  is  interactive  and  prompts  the 

operator  to  input  values  for  L,  M,  SNR,  DSEED,  and  a .  The 

operator  is  then  prompted  to  enter  values  for  the  pupil 

vector  and  object  vector.  The  program  uses  the  pupil  vector 

-k  — k 

to  generate  the  P*P  vector  which  in  turn  provides  the 
elements  of  the  diagonal  matrix,  B.  The  IMSL  subroutine 
VCONVO  performs  the  necessary  convolution.  The  values  for  M 
and  N  are  sufficient  to  create  the  matrices  F  (Fourier 
transform  matrix),  FF  (inverse  Fourier  transform  matrix),  and 
D  (spatial  truncation  matrix  of  M  diagonal  elements).  With 


the  matrix  B  available,  the  matrix  product  (FF) (B) (F) (D)  =  A 
(T  in  program  SRES)  can  be  calculated.  After  the  matrix  A  is 
calculated,  the  product  Axo  can  be  calculated  to  yield  the 
noise  free  image.  With  the  noise  free  image  available,  the 
noisy  image  can  be  calculated. 

The  operator  has  already  input  the  value  for  the  signal 
to  noise  ratio  (SNR)  which  can  be  used  with  the  noise  free 
image  to  create  a  noisy  image.  The  noise  is  assumed  to  be 
Gaussian,  white,  and  additive  for  all  spatial  frequencies. 

The  definition  of  SNR  is 

SNR  =  PSEQ  (4.1) 

On 

which  means  that 

On  =  PSEQ  (4.2) 

SNR 

where  the  noise  variance  is  represented  as  On  and  PSEQ  is  the 

— * 

power  in  the  image  sequence.  The  power  in  the  vector  x  is 
defined  as 


M 

PSEQ  =  (1/M)  £  (xk)J  (4.3) 

k=l 

where  the  vector  is  of  length  M.  The  noise  variance,  o„  •  is 
the  power  in  the  noise  vector  xn  and  is  defined  as  shown  in 
equation  (4.3).  The  noise  variance  assumes  a  zero  mean.  The 
variance  is  also  equal  to  the  square  of  the  standard 


deviation  of  the  noise  distribution.  With  SNR  being  provided 
by  the  operator,  the  variance  can  be  calculated  from  equation 
(4.2).  Using  the  IMSL  subroutine  GGNML  and  the  value  of  the 
variance,  the  Gaussian  noise  vector  is  generated.  The  noise 
is  amplitude  at  this  point,  so  each  term  is  squared  for  power 
and  added  termwise  to  the  original  noiseless  image  to  create 
the  noisy  image.  At  this  point,  the  least  squares  process 
can  be  implemented.  After  the  matrix  A  (T  in  computer 
program)  has  been  created,  the  IMSL  subroutine  LSVDF  is  used 
to  create  the  three  matrices,  CJ,S,V',  which  are  to  be  used  in 
the  superresolution  process.  After  the  singular  values  have 
been  generated,  the  substitution  for  l/sk  is  implemented 
using  the  smoothing  parameter,  » ,  to  negate  the  effects  of 
the  very  small  singular  values  of  A. 

As  can  be  easily  seen,  the  computer  program  performs  the 
exact  operations  described  in  chapter  III.  In  the  next 
chapter,  the  results  for  various  pupil  configurations  and 
parameter  values  will  be  examined.  It  will  be  shown  that  the 
least  squares  method  with  a  priori  knowledge  can  resolve  two 
incoherent  point  sources  not  ordinarily  resolved  in  a 
conventional  low  pass  imaging  system. 
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Figure  5.  Superresolution  Flow  Chart 


V.  Results  and  Conclusions 

Introduction 

This  chapter  presents  the  results  and  conclusions 
obtained  using  the  computer  program  SRES  to  simulate 
superresolution  using  the  least  squares  method  outlined  in 
the  previous  chapters.  In  review,  1/2  P  is  one-half  of  the 
symmetric  pupil  function  of  length  L,  where  L  <=  16.  M  is 
the  assumed  or  known  maximum  length  of  the  object,  where  M  <= 
32  since  N=32  is  the  square  dimension  of  the  transfer  matrix 
A  (T  in  SRES) .  SNR  is  the  signal  to  noise  ratio  and  is 
the  smoothing  or  regularization  constant  used  to  attenuate 
the  effect  of  the  very  small  singular  values.  The  error  used 
throughout  this  chapter  is  the  Euclidean  distance  or  two  norm 
of  the  difference  vector  between  the  object  and  noisy  image 
where  the  image  is  the  superresolved  image.  The  effects  of 
the  different  parameters  will  be  investigated  and  discussed. 
After  the  results  are  presented,  a  conclusion  section  will 
close  out  the  formal  portion  of  this  thesis.  Also,  pupil 
function  performance  for  various  pupils  is  contained  in 
appendix  B. 

Results 

Band-Pass  vs  Low-Pass  Pupil .  Both  pupils  considered 
contained  six  apertures  corresponding  to  the  six  elements  of 
the  pupil  vector.  The  low  pass  half-pupil  was  (OOOOOIII)', 
and  the  high  pass  half-pupil  was  (11100000)'.  Both  pupils 
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column  vectors  as  denoted  by  the  transpose  symbol .  The 
modulation  transfer  functions  ( MTF )  for  each  of  these  pupils 
is  shown  in  Fig.  6.  As  can  be  seen,  the  band-pass  pupil 
passes  higher  frequency  components  than  the  low-pass  pupil. 
The  band-pass  pupil  passes  the  higher  spatial  frequency  terms 
needed  to  resolve  the  two  incoherent  point  sources  located 
closer  than  the  normal  resolution  distance.  Figure  7 
illustrates  the  effect  of  using  the  high-pass  pupil  instead 
of  the  low-pass  pupil  for  superresolution.  By  using  the 
band-pass  pupil,  the  superresolution  scheme  served  to 
separate  the  two  incoherent  point  sources,  with  only  two 
small  (40%)  side  lobes  present.  By  using  a  threshold 
detection  criteria,  the  high  pass  pupil  can  superresolve  the 
object  consisting  of  the  two  incoherent  point  sources  using 
the  least  squares  method  and  a  smoothing  function,  while  the 
low  pass  pupil  cannot  superresolve  the  object. 

Effect  of  A  Priori  Knowledge .  Figure  8  illustrates 
the  importance  of  a  priori  knowledge  using  the 
superresolution  algorithm.  By  using  the  low-pass  pupil 
shown,  the  normal  image  does  not  resolve  the  two  incoherent 
point  sources  (object)  of  Figure  7-a) .  The  superresolution 
scheme  is  able  to  resolve  the  object  only  at  M=8,  where  M  is 
the  a  priori  knowledge  that  the  object  was  less  than  or  equal 
to  8.  The  effects  are  quite  interesting  for  M=16,  as  the 
algorithm  tries  to  restore  the  object,  but  there  is 
insufficient  information  to  do  so.  Figure  11  presents  a  plot 
of  M  vs  Error  for  the  five  various  pupil  functions  at  a 


35 


insufficient  information  to  do  so.  Figure  11  presents  a  plot 
of  M  vs  Error  for  the  five  various  pupil  functions  at  a 
constant  SNR  =100.  As  can  be  seen,  the  known  information 
concerning  the  object  does  reduce  the  error  in  many  cases. 

It  is  interesting  to  notice  that  in  some  cases,  the  error  is 
almost  constant  for  the  values  of  M  up  to  the  bend  in  the 
curve.  The  error  for  P3,  P4,  and  P5  are  all  fairly  constant 
up  to  M=24 .  The  smaller  pupils  exhibit  similiar  responses 
for  smaller  values  of  M.  It  is  important  to  define  the 
object  domain  as  close  as  possible  in  order  to  achieve 
reliable  and  accurate  results. 

Effect  of  a.  Figure  9  illustrates  the  effect  of  the 
smoothing  parameter,®.  As  can  be  seen,  holding  all 
parameters  but  constant,  the  higher  values  for 
( 10E-5 , 10E-10)  serve  to  resolve  the  object  quite  well  with  a 
rather  limited  pupil.  As  a aproaches  zero,  the  effect  of  the 
noise  due  to  the  ill-conditioned  nature  of  the  system  is 
obvious.  As  can  be  seen,  the  processed  images  using  very 
small  values  for  ®  are  quite  haphazard  and  look  nothing  like 
the  original  object,  even  with  the  high  SNR=100.  Due  to  the 
nature  of  the  imaging  problem,  the  superresolution  alogorithm 
must  use  a  good  estimate  for  the  smoothing  parameter,  a  ,  in 
order  to  achieve  accurate  results. 

Effect  of  Noise .  Figure  10  serves  to  illustrate  the 
effects  of  noise  on  the  superresolution  technique  using  the 
least  squares  method  with  a  smoothing  parameter.  The 


superresolved  images  all  appear  to  be  the  same,  illustrating 
the  noise  resistant  nature  of  the  superresolution  algorithm. 
Figure  12  also  shows  the  noise  resistance  of  the  least 
squares  method.  For  SNR  equal  to  100,  50,  and  5;  the  error 
vs  pupil  size  is  almost  equal  for  each  pupil.  Only  at  the 
high  SNR  of  2  does  the  noise  dominate.  This  algorithm  has 
proved  to  be  quite  good  in  the  presence  of  moderate  noise. 

Error  Analysis .  Figures  10  and  11  illustrate  the  sources 
of  error  present  using  the  least  squares  method  for 
superresolution.  The  error  is  expressed  as  the  distance 
between  the  object  and  superresolved  image,  using  the 
Euclidean  distance,  or  two  norm. 

In  Figure  10,  the  error  is  plotted  on  the  vertical  axis 
and  the  pupils  are  plotted  on  the  horizontal  axis.  In  this 
graph,  all  pupils  are  the  high  pass  versions  with  PI 
consisting  of  6  elements,  P2  consisting  of  8  elements  and  so 
forth  with  the  last  pupil,  P5  consisting  of  14  elements. 

Each  curve  is  at  a  constant  SNR.  The  error  is  seen  to 
decrease  as  the  pupil  elements  are  increased.  This  is 
logical,  since  the  more  pupil  elements  available,  the  more 
spatial  frequencies  passed  by  the  system,  resulting  in  more 
information  available  for  the  algorithm. 

In  Figure  11,  the  effect  of  a  priori  information  is 
analyzed.  The  known  extent  of  the  object,  M,  is  plotted 
versus  error  for  constant  pupils.  At  the  extreme,  with  the 
object  known  to  be  within  32  bits,  the  larger  pupils  have  the 
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a)  Processsed  Image  M=32  b)  Processed  Image.  M=32 

«=  10E-10 ,  SNR=100  *  =  10E-10 ,  SNR=50 

1/2P  =[lllOOOOO]  1/2  P  =[11100000] 


c)  Processed  Image,  M=32 
or  =  10E-10,  SNR=5 . 0 
1/2  P  =[11100000] 


Figure  10.  Superresolution  Example  #4 


least  error,  which  Is  logical.  As  the  known  extent  of  the 
object  Increases,  or  as  M  decreases,  the  error  decreases, 
with  the  larger  pupils  approaching  the  constant  level  quite 
rapidly.  For  the  10,  12,  and  14  bit  pupils,  the  error  is 
constant  up  to  M=24,  while  the  6  and  8  bit  pupils  decrease  at 
different  rate  that  appears  almost  linear.  The  effect  of  M 
is  obvious,  but  as  can  be  seen  from  some  of  the  pupils,  the 
increase  in  information  concerning  the  maximum  extent  does 
not  yield  appreciable  results  concerning  error  reduction.  As 
an  example,  for  P3  consisting  of  10  bits,  the  error  for  M=24 
is  about  the  same  as  for  M=8 .  For  P3,  the  information  gained 
by  knowing  that  M=8  instead  of  at  least  24  is  slight.  The 
error  is  seen  to  be  almost  linear  over  some  range,  and  almost 
flat  up  to  a  threshold  level. 

Conclusions 

Based  on  the  performance  of  the  computer  simulation, 
superresolution  is  achievable  in  limited  bandwidth  optical 
systems.  The  performance  of  the  superresolution  process  is 
affected  by: 

A  priori  information  available  concerning  the  object 

Type  of  pupil 

And  noise  in  the  system. 

By  using  a  smoothing  parameter,  a  ,  the  overamplification  of 
the  high  order  noise  terms  is  avoided.  With  finite  precision 
arithmetic  and  noise,  the  very  small  singular  values  of  the 
transfer  matrix  must  be  attenuated  in  order  to  preserve  the 
identity  of  the  original  object.  The  problem  with 


attenuating  the  high  order  singular  values  is  that  these 
singular  values  correspond  to  the  high  order  spatial 
frequency  terras  needed  for  resolution.  Obviously  some  high 
order  terms  are  necessary  for  resolution,  but  amplification 
of  the  noise  breaks  down  the  superresolution  process  when 
high  order  singular  values  are  present.  For  this 
one-dimensional  case,  the  pupil  was  assumed  to  be  symmetric 
and  quite  limited  with  a  maximum  length  of  16.  For  a  more 
rigorous  approach,  the  pupil  length  could  be  increased,  but 
the  resulting  increase  in  computer  complexity  would  be  quite 
costly.  The  results  obtained  graphically  and  in  Appendix  B 
illustrate  the  importance  of  the  different  parameters  and 
verify  the  premise  that  high  pass  or  band  pass  pupils  will 
image  better  with  this  particular  superresolution  algorithm 
than  low  pass  pupils  of  the  same  bandwidth.  By  knowing  as 
much  as  possible  about  the  noise  and  object,  superresolution 
is  achievable  and  possible  using  the  least  squares  method. 
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C  REPLACE  SCK)  WITH  F(S(K)) 

C  THIS  STEP  REPLACES  SCK)  WITH  FCS(K)) 

C  ANC  PERFORMS  THE  COT  PRCCLCT 

C  WITH  THE  LEFT  SINGULAR  VECTCfi  FRCM  CCLS  CF  LT 

C  SAVE  SINGLLAR  VALUES  AS  SV 

CO  136  K*  1  »  M 

SV(K)=SCK  > 

>  X(K)  =  O.C 
CC  134  I  5 1  ♦  N 

>X(K)=XCC(I)*UTCK,I)4XXCK) 

134  CCNTINLE 

SS(K)*CS(K)*S<KXS(K))/(CS(lO*SCK)*S(K)*SCK))-*ALFHA) 

SS(K)*SS(K)*XX(lO 

136  CONTINUE 

C  MULTIPLY  COT  PRCCLCT  X  RIGHT  SINGULAR  VECTORS  FRCM  T 

C  YYCI)  IS  RECONSTRUCTED  IMAGE 

DO  137  I s 1 » M 

XXCIXO.C 
CC  138  J  *  1  *  M 

XX(I)*XXCI)4TCI,J)*SSCJ) 

138  CONTINUE 

YY(I)*XX(I) 

137  CONTINUE 

C  NORMALIZE  RECONSTRUCTED  IMAGE 

CC  162  1*1, M 

IFCYY(I).NE.C.O)S(l)=ABSCYYCI)) 

IF  CSCD.GT.C.O)  GCTC  163 

162  CONTINUE 

163  CC  164  1*1,  M 
IF(A6S(YY(I)).GT.S(1))S(1)sAES(YY(I)) 

164  CCNTINLE 

C  SCI)  IS  MAX  AES  VALUE  OF  R  EC  IMAGE 

CO  167  1*1, M 

YYCI)*YYCI)/S(1) 

167  CCNTINLE 

WRITE  <96  *1775 
W»ITE  (96,178) 

WRITEC96.20) 

CC  139  1*1, N 
IF(I.G7.M)YY(I)*C.O 
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IF(I.GT.LKCCI>  =  C.O 
IF(I.GT.N)SV(I)*C.O 
IF  (I.GT.LAUB)A(I>«C.O 

kaiTE(96,142>I,SVCI)«CwCl'),A<I),XCCI),XI(I),RK(I),X0C<n,YYCI) 

COMIME 

SS(1)=0.C 

XX(1)=0.C 

CC  146  1=1, N 

XX(1)=XX(1)*(YY(1)-XC(I))**2 
S$(1)=SS(1)  ♦  CXCC<X)-XCCX))*42 
CONTINUE 

SSC1)  IS  THE  EUCLIOEAN  DISTANCE  BETWEEN  NCISY 
CSSERVEC  IKAGE  XCC(I)ANC  NOISE  FREE  OBJECT ,  XCCI) 

XXC1)  IS  THE  EUCLICEAN  DISTANCE  BETWEEN  THE 
RcSTCREC  IHAGE  YY(I)ANC  THE  ORIGINAL  OBJECT  XCCI) 
SSC1)=SCRTCSSC1)) 

XX(1)=SCRT(XX(1)) 

WRITE  (9c, 20) 

V,  R I T  E  C  9  6  ,  2  C  ) 

WRITE  (96,209)SS<1) 

WRITE  (96,149)XX(1) 

WRITE  (96 , 17)SNR 
WRITE  (96 ,4)ALPHA 
WRITE(96,25)CCSEEC 


Appendix  B:  Data  from  Computer  Program  SRES 


Notes : 

1.  All  objects  are  two  incoherent  point  sources  two 
bits  apart 

2.  Processed  image  error  (PIE)  is  the  Euclidean 
distance  (two  norm  of  the  difference)  between  the 
superresolved  image  and  noise  free  object. 
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Appendix  C:  Typical  Data  from  SRES  computer  Program 


SV  =  Singular  Value 

P*P=  Autocorrelation  of  Pupil  Vector 
New  Image  is  processed  image 

Error  is  Euclidean  distance  between  the  image  and  object 
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